In [1]:
%matplotlib inline
%config InlineBackend.figure_format = "retina"
from matplotlib import rcParams
rcParams["savefig.dpi"] = 100
rcParams["font.size"] = 20
In [2]:
import fitsio
import numpy as np
import matplotlib.pyplot as pl
In [3]:
times = fitsio.read("centroids.fits")
data = fitsio.read("centroids.fits", 1)
Remove the stars that have NaNs at all times.
In [280]:
m = np.sum(~np.any(np.isnan(data), axis=1), axis=-1) > 10
for k in range(len(data)):
if m[k]:
m[k] *= np.all(np.var(data[k][:, np.all(np.isfinite(data[k]), axis=0)], axis=1) < 1.0)
good_stars = data[m, :, :]
In [281]:
print(good_stars.shape)
Remove times that have NaNs in all pixels.
In [282]:
m = np.all(np.isnan(good_stars), axis=(0, 1))
X = good_stars[:, :, ~m]
good_times = times[~m]
In [283]:
mask = ~np.any(np.isnan(X), axis=1)
K, T = mask.shape
In [401]:
nlatent = 3
tt = (good_times - np.mean(good_times)) / np.std(good_times)
W = np.vander(tt, nlatent + 1)[:, ::-1]
W[:, 1:] = (W[:, 1:] - np.mean(W[:, 1:], axis=0)) / np.std(W[:, 1:], axis=0)
A = np.empty((2 * K, nlatent + 1))
Update all the distortions given random pointings.
In [409]:
for k in xrange(K):
A[2*k:2*k+2, :] = 2
Use this A to update W.
In [408]:
for t in xrange(T):
mt = (mask[:, t][:, None] + np.zeros((K, 2), dtype=bool)).flatten()
xt = X[:, :, t].flatten()[mt] - A[mt, 0]
W[t, 1:] = np.linalg.solve(np.dot(A[mt, 1:].T, A[mt, 1:]) + 1.0 * np.eye(nlatent), np.dot(A[mt, 1:].T, xt))
W[:, 1:] -= np.mean(W[:, 1:], axis=0)
u, s, v = np.linalg.svd(W[:, 1:])
W[:, 1:] = np.dot(v, W[:, 1:].T).T
In [410]:
np.mean([np.mean((X[k, :, mask[k]] - np.dot(A[2*k:2*k+2], W[mask[k]].T).T) ** 2) for k in xrange(K)])
Out[410]:
In [416]:
pl.figure(figsize=(10, 10))
for i in range(nlatent):
pl.subplot(nlatent, 1, i+1)
pl.plot(good_times, W[:, i+1], ".k")
pl.xlim(1960, 1970)
In [415]:
pl.plot(good_times, X[0, 0, :], "b")
pl.plot(good_times, X[0, 0, :], ".k", ms=2)
pl.xlim(1960, 1970);
In [418]:
fits = fitsio.FITS("pointing_model.fits", "rw")
fits.write(good_times)
fits.write(W)
fits.close()
In [ ]: